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Abstract 

We describe rabacus, a Python package for calculating the transfer of hydrogen ionizing radiation in simplified geometries relevant 
to astronomy and cosmology. We present example solutions for three specific cases: 1) a semi-infinite slab gas distribution in a 
homogeneous isotropic background, 2) a spherically symmetric gas distribution with a point source at the center, and 3) a spherically 
symmetric gas distribution in a homogeneous isotropic background. All problems can accommodate arbitrary spectra and density 
profiles as input. The solutions include a treatment of both hydrogen and helium, a self-consistent calculation of equilibrium 
temperatures, and the transfer of recombination radiation. The core routines are written in Fortran 90 and then wrapped in Python 
leading to execution speeds thousands of times faster than equivalent routines written in pure Python. In addition, all variables have 
associated units for ease of analysis. The software is part of the Python Package Index and the source code is available on Bitbucket 
at https: //bitbucket. org/galtay/rabacus In addition, installation instructions and a detailed users guide are available at 
http://pyth 0 nh 0 sted. 0 rg//rabacus 
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1. Introduction 

The transport of radiation is central to model building in 
astronomy and cosmology, yet few closed form solutions ex¬ 
ist even in highly symmetric situations. In this work we de¬ 
scribe rabacus, a Python package for calculating equilibrium 
ionization and temperature states for well posed radiative trans¬ 
fer problems in planar or spherical geometries. 

There is a long history of similar numerical approaches in 
the field of photon dominated regions or photo-dissociation re¬ 
gions (PDRs, see Hollenbach and Tielensj [1999 Rollig et al 


2007 for reviews). PDRs are usually defined as regions of the 
interstellar medium in which far ultraviolet photons with ener¬ 
gies between 6 and 13.6 eV dominate the thermal and chemical 
dynamics. Typical densities and temperatures are n u > 10 cm -3 
and T < 5000 K. As the main goal of PDR codes is the calcula¬ 
tion of emission spectra, they typically include large chemical 
networks and take account of many line excitation processes. 

The goal of rabacus is to provide a similar tool that is op¬ 
timized for the transfer of hydrogen ionizing radiation through 
cosmological scale systems. Common applications include quasar 
absorption systems ( Wolfe et al.| 2005 1 Meiksin 2009|), the 


evaporation of mini-halos during reionization (Shapiro et al. 


2004), and Stromgren spheres ( [Stromgren] 1 1939| ). Typical den¬ 
sities and temperatures in these regions are n H < 10 _1 cm -3 
and T > 8000 K. This allows us to consider a compact chem¬ 
ical network consisting of hydrogen and helium and limits the 
list of physical processes that dominate the dynamics to photo¬ 
ionization, photo-heating, collisional ionization, collisional ex¬ 
citation, recombination, free-free radiation, Compton scatter¬ 
ing, and Hubble expansion. 

The density and temperature regime treated by rabacus es¬ 


sentially excludes the cold inter-stellar medium phase of gas 
where molecules and dust begin to play an important role in the 
thermal and chemical state of hydrogen and helium. However, 
cooling from metal line emission can play an important role 
in diffuse enriched gas (e.g. Wiersma et al.||2QQ9 ). Metal line 
cooling is not treated by rabacus. However, its effect on hy¬ 
drogen and helium can be approximated by first calculating the 
equilibrium temperature of a system without metals and then 
running iso-thermal models with lower temperatures. 

Photo-ionization and heating rates are calculated in rabacus 
using a multi frequency reverse ray tracing technique ( |Altay| 
|and Themis] |2013| ). In cases of planar geometry this is sim¬ 
ply a matter of evaluating exponential integrals, but in the case 
of spherical geometries embedded in isotropic radiation back¬ 
grounds, ray casting is used. This allows us to use the exact path 
length through each shell in our radiative transfer solutions. 

The development of rabacus has several purposes. One is 
the treatment of geometries and radiation fields that are not in¬ 
cluded in many popular codes in use today (for example plane 
parallel radiation incident from both sides of a slab or isotropic 
backgrounds incident on spherical geometries). In addition we 
wanted to use an object oriented approach coupled with a mod¬ 
ern language with a large user base and fully developed graphi¬ 
cal libraries. In addition, we wanted to make the isolation of 
physical processes convenient. We give special attention to 
physical processes that are most often treated in an approximate 
way in cosmological hydrodynamic simulations: 1) the spec¬ 
tral shape of the ionizing background, 2) the effect of helium 
on spectral hardening and hydrogen ionization profiles, 3) the 
effect of equilibrium temperatures vs. fixed temperatures, and 
4) the effect of using the on-the-spot approximation vs. calcu¬ 
lating the transfer of recombination radiation (e.g. Altay et al 
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201 It McQuin n et al.[ |2011[ |Friedrich et al.| |2012| |Rahmati| 
et al.[|201?||Raicevic et al.[|2014| ). 

We developed the Python package rabacus with the goal of 
making these solutions both available and readily usable by the 
astrophysics community at large. To this end, our project is 
free, open source, and makes use of a high level programming 
language with a large user base and well supported graphics 
packages, rabacus is available on the Python Package IndesQ 
the source code is available on Bitbucke0 and a users guide 
is available through pythonhosted.or^ The software is re¬ 
leased under the terms of the FreeBSD licens^] While Python 
has many benefits, the fact that it is an interpreted language 
can make computationally intensive tasks slow compared to 
compiled codes. In order to circumvent this problem we im¬ 
plemented the core functionality of rabacus in Fortran 90 and 
wrapped the resulting code in Python using the Numpy package 
f2py. In addition, all variables in rabacus have units through use 
of the Quantities Python package^] This package will automat¬ 
ically be downloaded during installation. Comprehensive ex¬ 
amples in the users guide cover the use of variables with units. 
The rest of this paper is organized as follows. In §2 we present 
our notation and review some basic physics. In §3 we discuss 
classes that represent sources of radiation. In §4 we describe 
the single zone solvers supplied by rabacus. In §5-8 we discuss 
the geometric solvers, and in §9 we conclude. 


2. Basic Physics and Notation 


In this section we review the physics governing the temper¬ 
ature and ionization state of cosmological gas. We will consider 
the number density of six species and indicate ionization states 
using Roman numerals { n m , n mi , n HeI , n Hen , n Hem , n e }. The free 
electron number density is related to the other species as, 
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and the mass fraction of helium is given as, 
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In what follows, we make use of the following notation, 
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where the electron fraction, v , is the ratio of the number den- 

7 e 1 

sity of electrons and the maximum possible number density of 
electrons, 
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1 https://pypi.python.org/pypi/rabacus 

2 https://bitbucket.org/galtay/rabacus 

“http://pythonhosted.org/rabacus 

4 http://opensource.org/licenses/BSD-2-Clause 

5 https://github.com/python-quantities/python-quantities 
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2.7. Ionization Evolution 

Considering the processes of photo-ionization, collisional 
ionization, and recombination, the coupled evolution equations 
for the ionization fractions take the following form, 
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where the T’s represent photo-ionization rates, the C’s col¬ 
lisional ionization rates, and the R’s recombination rates. Col¬ 
lisional ionization and recombination rates have a non-trivial 
dependence on temperature and rabacus provides a simple in¬ 
terface to load the rate fits presented in |Hui and Gnedin) \\991) . 
Photo-ionization rates depend on the spectrum of ionizing radi¬ 
ation. 
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where u v is the energy density per unit frequency in radiation, r y 
is the frequency dependent optical depth in all species, the N x 
are column densities, the <x x are frequency dependent photo¬ 
ionization cross sections, and the v* are ionization thresholds. 
Note that u v represents an optically thin quantity and exp(-r y ) 
accounts for absorbing material between the source of radiation 
and the point at which the photo-ionization rate is being calcu¬ 
lated. 

rabacus provides tools for the creation of spectra and access 
to the cross section fits of |Vemer et al.||1996| ). The variable u v 
serves to characterize the radiation field in a geometry indepen¬ 
dent way. It is related to more geometry specific characteriza¬ 
tions as follows, 


CUy = 4 71 Jy = Fy = (9) 

4/rr z 

where J v is the angle averaged specific intensity, F v is the en¬ 
ergy flux per unit frequency of a plane parallel source, L v is the 
energy luminosity per unit frequency of a point source and r is 
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the distance from the point source. It is sometimes convenient 
to group terms using the following notation, 


Their contribution to the cooling function is calculated as fol¬ 
lows, 
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However, it is important to keep in mind that n G depends on the 
ionization fractions. Using this notation we can express these 
equations in matrix form, 
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2.2. Temperature Evolution 

The evolution of gas internal energy u is governed by the 
following equation, 


3 

« = ^k b T(n H +n He +n e ) (13) 

du 

— - <H- A c -3Hk h T(n u + n Ue + n e ) (14) 


The first term, *H, accounts for photo-heating, the second term, 
A c , accounts for radiative cooling, and the final term is due to 
adiabatic cooling from Hubble expansion. The first two are al¬ 
ways included in the rabacus solvers, but the last is optional. 
The photo-heating term depends on the spectrum of the radia¬ 
tion field. 
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The first nine terms are products of temperature dependent rates 
and ion number densities. The next term accounts for bremsstrahlung 
cooling and involves a Gaunt factor which we set to gff = 1.1 + 

0.34 exp[-(5.5-log 10 T) 2/3 ]. The final term accounts for Comp¬ 
ton cooling. The fits used in rabacus for these rates are from 
|Hui and Gnedin] \\991) . 


3. Source Classes 

There are three source classes available in rabacus namely, 
PointSource, PlaneSource, and BackgroundSource. The 
sources are grouped by geometry due to the different kinds of 
integrals that need to be computed in order to calculate quan¬ 
tities like photo-ionization and photo-heating rates. All classes 
presented here take a common set of arguments when instan¬ 
tiated. The first two, qjnin and qjnax are floats equal to the 
minimum and maximum photon energies considered divided by 
Rydbergs. The third determines the shape of the spectrum and 
is a string chosen from the following list { monochromatic, 
hml2, thermal, powerlaw, user } where hml2 is the spectral 
model from |Haardt and Madau| ( |2012] ) (HM12) and user is a 
user defined spectrum. Keyword arguments supply extra infor¬ 
mation if needed. For example, the effective temperature in the 
case of thermal spectra, the powerlaw index in the case of pow¬ 
erlaw spectra, or the redshift in the case of hml2. All source 
classes contain geometry specific normalization routines. For 
example, a point source can be normalized to produce a given 
luminosity in erg s -1 while a plane source can be normalized to 
produce a given flux in erg cm -2 s -1 . 

By default, if the range of energies requested spans the hy¬ 
drogen and helium ionizing thresholds (i.e. if qjnin < 1 and 
q_max > v* /v* ~ 4) then frequency averaged, or “grey”, 
photo-ionization cross sections will be calculated as follows, 


We consider the following radiative cooling processes, 

• Collisional Ionization Cooling (cic) 

• Collisional Excitation Cooling (cec) 

• Dielectronic Recombination Cooling (di) 

• Radiative Recombination Cooling (re) 

• Bremsstrahlung Cooling (br) 

• Compton Heating/Cooling (cp) 
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where v max = q_max x vjjj, and qjnax is one of the arguments 
provided to the source class. This allows for implicit definitions 
of grey frequencies and energies, 

cr x (vf y ) = crf y , £f y = /zvf y (18) 


where X is one of { HI, Hel, Hell}. However, these quantities 
are calculated only as a convenience to the user. Unless a mono¬ 
chromatic spectrum is specified, the geometric solvers provided 
by rabacus, described in perform multi-frequency radiative 
transfer. 

By default, all spectra are sampled at N v = 128 photon en¬ 
ergies logarithmically spaced between q_min and q_max. The 
variable N v can be passed into the source initialization rou¬ 
tines to increase or decrease spectral resolution. All photo¬ 
ionization and photo-heating rates are calculated by performing 
discretized versions of the integrals in Eqs. [8] and |T5] Photo¬ 
ionization cross-sections are calculated by evaluating the fitting 
functions of Verner et al.| (1996| at the N v photon energies. User 
defined spectra are defined by providing the energy samples 
and spectral shape in tabulated form. In general, increasing the 
value of N v increases the accuracy of the integrations at the cost 
of longer computations. 


4. Single Zone Solvers 

The single zone solvers provide equilibrium solutions to 
Eqs. m and [14] Three classes are provided by rabacus to 
handle three types of equilibrium. 


4.1. Collisional Equilibrium 

The first case we consider is that of collisional equilibrium 
handled by the class Solve_CE. In this case, temperatures are 
fixed, photo-ionization rates are zero, and ionization fractions 
are found such that collisional ionizations balance recombina¬ 
tions. In this case, there are closed form solutions for the ion¬ 
ization fractions, 
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These solutions also determine the minimum possible electron 
density for a given temperature, n™ m = x^n n + (x^ n +2 )n He . 


4.2. Photo Collisional Equilibrium 

The second case we consider is that of photo-collisional 
equilibrium handled by the class Solve_PCE. In this case, tem¬ 
peratures are fixed, photo-ionization rates are finite, and ioniza¬ 
tion fractions are found such that collisional and photo ioniza¬ 
tions balance recombinations. In this case, implicit solutions 
for the ionization fractions are, 
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They are implicit solutions due to the factors of n Q on the 
right hand side. It is possible to eliminate the n e yielding a 
univariate fourth order polynomial equation for each ionization 
fraction, however this yields solutions with thousands of terms. 
In addition, the floating point errors involved in numerically 
evaluating these solutions with fixed precision frequently pro¬ 
duce imaginary components. 

In rabacus we use a robust iterative approach. It involves an 
initial guess for n e using a closed form solution for the hydrogen 
ionization fraction (see §6.1). This allows for the evaluation 
of all ionization fractions using Eqs. ]2l] and [22] These are 
then used to re-calculate n t using Eq. gjwhich in turn can be 
used to re-calculate the ionization fractions. This process is 
repeated until the amplitude of changes in n e drops below a 
given tolerance. In most cases, the number of floating point 
operations for the iterative method will be less than the number 
required to evaluate the thousands of terms in the closed form 
solution. In addition, the initial guess for n Q is bounded between 
a maximum value of n™ x = n u + 2 n Ue and a minimum value 
given by collisional ionization equilibrium. 


4.3. Photo Collisional Thermal Equilibrium 

The third case we consider is that of photo-collisional-thermal 
equilibrium handled by the class Solve_PCTE. In this case, ion¬ 
ization fractions are found such that collisional and photo ion¬ 
izations balance recombinations, and a temperature is found 
such that cooling balances heating. As described above, we use 
an iterative technique to find solutions. First a temperature is 
guessed then photo-collisional equilibrium is found at that tem¬ 
perature. Next, heating and cooling rates are calculated using 
the ionization fractions. If the heating rate is greater than the 
cooling rate, the temperature is increased and vice versa. This 
procedure is iterated until the heating rate is within a given tol¬ 
erance of the cooling rate. The range of temperatures to search 
can be defined by the user but defaults to T mm = 7.5 x 10 3 K 
and r max = 10 5 K. The initial temperature guess is logT = 
(log Unin + log T max )/2. 
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5. Geometric Solver Arguments 


In the next two sections we present the geometric solver 
classes available in rabacus and compare our solutions to those 
in the literature. These classes take a common set of arguments 
when instantiated. The first, Edges, defines the gas geometry. 
We focus on two geometries: spheres and semi-infinite slabs. 
Spherical gas distributions are centered on the origin and dis¬ 
cretized into spherical layers separated by infinitesimally thin 
spherical shells. In this case, Edges is an array of radial dis¬ 
tances indicating the position of each shell and typically has 
zero as the first element and the radius of the sphere as the last 
element. Slab gas distributions are discretized into rectilinear 
layers separated by infinitesimally thin planes. In this case, 
Edges is an array of distances indicating the position of each 
plane relative to the surface of the slab and typically has zero 
as the first element and the thickness of the slab as the last ele¬ 
ment. For N discrete layers, Edges will have N+l entries. The 
next three arguments, T, nH, nHe, are arrays of length N describ¬ 
ing the temperature, number density of hydrogen and number 
density of helium in each layer. 
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where z is depth into the slab, cr is the hydrogen photo-ionization 
cross-section at the monochromatic frequency considered, is 
the neutral fraction at the surface of the slab, and x™ is the neu¬ 
tral fraction in collisional equilibrium at the fixed temperature 
of the slab. Note that this is an inverse solution in the sense that 
it provides depth into the slab as a function of neutral fraction. 
However, it allows the full ionization profile to be mapped out 
by plugging in values of x m between the two extremes x^ and 
is a closed form solution for x ce . There is also a 


Eq. 19 
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The coefficients in this quadratic solution are, 


6. Slab Solvers 
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rabacus has three slab classes: SlabPln which models plane 
parallel radiation incident from one side; Slab2Pln which mod- 
els plane parallel radiation incident from both sides; and Slab2Bgnd where y represents the contribution to free electrons from el- 


which models a slab in a uniform isotropic background. Our 
first goal is to validate the rabacus slab solver using a closed 
form solution in the case of monochromatic radiation and uni¬ 
form density and temperature. Afterward, we demonstrate the 
effects of polychromatic spectra, temperature evolution, and re¬ 
combination radiation. 

To begin, we create two plane parallel radiation sources us¬ 
ing the PlaneSource class. The first has a spectrum taken 


from the cosmological UV background model of Haardt and 
Madau|p012| ) (HM12) at z = 3 between 1 and 400 Rydbergs. 


We then use the neutral hydrogen grey energy (see Eq. [18]) to 
create a monochromatic instance of PlaneSource and normal¬ 
ize it such that it has the same hydrogen photo-ionization rate 
as the polychromatic instance of PlaneSource. This results 
in a source with photons of energy 16.8 eV, a photon flux of 
2.31 x 10 5 cm -2 s -1 , and a hydrogen photo-ionization rate of 
T hi = 8.27 x 10 -13 s -1 . We then follow the same procedure us¬ 
ing the BackgroundSource class. This produces four source 
objects (a pair of monochromatic sources and a pair of poly¬ 
chromatic sources used to find the grey energies) we can pass 
to the geometric solvers. 


ements other than hydrogen, n e = (x HU + y)n u , and T^ is the 


photo-ionization rate at the surface of the slab (Altay and The- 
|uns[ |2013] ). 

To generate an appropriate slab model in rabacus, we use 
the Cosmology class to calculate the critical density of hydro¬ 
gen at z = 3 using the Planck cosmological parameters ( |Planck| 
[Coll aboratio n et al.| |2013| ) and then calculate the Jeans length 
(assuming a mean molecular weight p - 1) for gas at 10 4 K 
and A = 2200 times the critical density. The value for A was 
arbitrarily chosen to produce a small neutral region in the slab 
and results in a length of L = 4.17 kpc and number densities of 
n H = 8.61 x 10 -3 cm -3 and n Ue = 7.09 x 10 -4 cm -3 . 

In Fig. |T] we show the closed form (x a ” a , x aaa ) and raba¬ 
cus solutions. We generated the rabacus solution by passing 
the monochromatic plane source and the slab description to the 
SlabPln class. The maximum relative and absolute difference 
between the two solutions occurs near the ionization front but is 
always less than 7 x 10 -3 and two orders of magnitude smaller 
than that over most of the length of the slab. This indicates that 
the basic machinery of the rabacus solver is correct. 

6.2. Isotropic Background, Recombination Radiation 


6.1. Plane Parallel Closed Form Solution 

In the case of monochromatic plane parallel radiation inci¬ 
dent from one side onto a pure hydrogen slab with uniform den¬ 
sity and temperature, there is a closed form solution describing 
the ionization profile if the on-the-spot approximation is used 
(see the Appendix of Altay and Themis] |20 13), 


In the previous section both the gas geometry and the radi¬ 
ation transport were one dimensional. In this section we dis¬ 
cuss isotropic backgrounds and diffuse recombination radiation 
both of which introduce angular variables into the problem. The 
Slab2Bgnd class models an isotropic background incident onto 
a planar gas distribution. Here we describe the equations that 
are solved by rabacus to account for both the background ra¬ 
diation and radiation from recombinations. We also note that 
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Figure 1: Ionization fractions for the closed form solution test described in 
section [6~T] In the lower panel we show ionization fractions from rabacus (thick 
lines) and the closed form solution (thin lines). In the middle panel we show the 
absolute difference in ionization fractions x ana - rabacus. In the top panel we 
show the relative difference in ionization fractions ( x ana - rabacus ) / rabacus. 
In the bottom panel, the rabacus solution is shown in black while the analytic 
solution is shown in green. 
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Figure 2: The geometry used in the SlabBgnd class. 


recombination radiation can optionally be treated the same way 
in the SlabPln and Slab2Pln classes. 

The angle averaged specific intensity at any depth z (see Fig. 
[2]) into the slab is, 

4(z)=TJ" d(p duly&fl) (26) 

When considering a point at depth z, it is convenient to decom¬ 
pose the rays into those with negative z-components (I~) for 
which id = cos 6 > 0 and those with positive z-components (/+) 
for which ji < 0. In addition, we will decompose the radiation 
into a component from the background itself (src) and a com¬ 
ponent from radiative recombinations in the slab (rec). Using 


this decomposition we can write, 
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where A r a ^ = r(a) - r(b), r(z f ) is the optical depth between 
z = 0 and z-z! along a path perpendicular to the surface of the 
slab, j y is the emission coefficient for recombination radiation, 
and is the un-attenuated isotropic background. The angle 
averaged specific intensity can now be written as, 



(31) 


The angular integrals can be expressed in terms of the exponen¬ 
tial integrals E\ and E 2 , 
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Combining all terms we get, 
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We account for recombination radiation from hydrogen and he¬ 
lium with monochromatic emission at the appropriate ioniza¬ 
tion threshold and calculate the emission coefficients as, 


Jhii 


7HeII = 


R „e„ H He A 'He„«e 


47T 


Ku n u X m r 


\n 


7HeIII 
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(37) 
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where the R l x represent recombinations to the ground state of 
each ion. In practice, this rate is calculated by subtracting the 
case B fit from the case A fit. This neglects some helium recom¬ 
binations that produce ionizing photons (see for example §3.1 
in Altay et al. 2008) but makes for a more straightforward com¬ 


parison to models which use case B rates only. Our treatment is 
similar to the models used in Faucher-Giguere et al. ( 2009| ) and 


Haardt and Madau ([ 2012 ) except that we include the effects of 


collisional ionization. 


6.3. Isolating Physical Processes 

Having verified the basic functionality of rabacus against 
a closed form solution, we now show the results of including 
more physical processes. All geometric solvers in rabacus in¬ 
clude the ability to treat polychromatic spectra, helium ioniza¬ 
tion, thermal balance, and the transfer of recombination radi¬ 
ation. We now use the slab geometry to produce five models 
which add these physical processes one at a time. For all of 
these examples we will use the Slab2Bgnd class which models 
a slab in an isotropic background of radiation. In general, plane 
parallel radiation will penetrate further into a slab than isotropic 
radiation. Therefore, we use a lower overdensity of A = 1800 in 
this test. This leads to a slab of length L = 4.61 kpc and number 
densities n H = 7.04 x 10 -3 cm -3 and n Ue = 5.80 x 10 -4 cm -3 . 
Although the temperature, and hence the Jeans length, will de¬ 
viate from 10 4 K when we solve for thermal balance, we keep 
the slab size fixed for these tests. This results in a neutral hy¬ 
drogen column density of N m = 5.67 x 10 18 cm -2 for our most 
realistic model. These volume and column densities are typi¬ 
cal for Lyman limit absorption systems. The five models we 
examine are, 

grey JixT_B monochromatic spectrum, fixed temperature, treats 
recombination photons by using case B rates 

hml2_fixT_B adds full HM12 spectrum as opposed to the mono¬ 
chromatic grey approximation 

hml2_evoT_B adds thermal balance as opposed to fixed tem¬ 
peratures 

hml2_evoT_A same as hml2_evoT_B but with case A rates 

hml2_evoT_rt uses case A recombination rates but adds the 
transfer of recombination photons from hydrogen and he¬ 
lium emitted isotropically from each layer 


the slab, N m , by a factor of 1.8. The monochromatic spectrum 
does not contain any helium ionizing photons and so a com¬ 
parison between the helium ionization profiles is not well moti¬ 
vated. 

By comparing models hml2_fixT_B and hml2_evoT_B we 
can isolate the effect of self-consistently calculating the tem¬ 
perature structure. The equilibrium temperature is greater than 
10 4 K at all depths in the slab and peaks at 1.6 xlO 4 K at the 
edge. This produces both higher collisional ionization rates and 
lower recombination rates leading to fewer neutral hydrogen 
and helium atoms. This effect is compounded by the additional 
photons that can penetrate into the more ionized slab. Calculat¬ 
ing the equilibrium temperatures lowers N m by a factor of 11 
and A HeI by a factor of 6.3 while A HeII is only affected at the five 
percent level. 

By comparing models hml2_evoT_B and hml2_evoT_A 
we get a sense of the maximum variation possible due to re¬ 
combination photons. In the case B model, recombination pho¬ 
tons are emitted and absorbed in the same layer, while in the 
case A model all recombination photons escape the system. In 
model hml2_evoT_rt we calculate the transfer of recombina¬ 
tion photons isotropically emitted from each layer. The rt and 
case A solutions agree at the optically thin surface of the slab. 
As depth increases, the rt solution diverges from the case A 
solution but always remains between the case A and case B so¬ 
lutions. This translates into a factor of 3.2 increase in both N m 
and A HeI relative to case B. These changes are larger than those 
that result from using a grey spectrum as opposed to a polychro¬ 
matic one. The columns of singly and doubly ionized helium 
are only affected at the few percent level. 

The fact that the rt solution does not converge to the case 
B result at any depth distinguishes the ionized slab model from 
models like the Stromgren sphere (see Q. The major differ¬ 
ence is that Stromgren spheres necessarily have a neutral re¬ 
gion enclosing the source of photons. Conservation of photons 
requires that all solutions agree at radii large enough where all 
photons have been absorbed. However, when the slab is not 
fully neutral at its center (as will be the case for Lyman limit 
systems as well as lower column density systems), the rt solu¬ 
tion will differ from the case B solution at all points in the slab. 
In summary, the behavior of the most realistic solution here is 
not captured by either case A or case B approximations. In ad¬ 
dition, the total column density in neutral hydrogen and helium 
differs between the models by at least a factor of 3 and at most 
a factor of 11. 


In Fig. [3] we show ionization and temperature profiles for all 
five models. Note that we have taken advantage of the symmet¬ 
ric geometry and only show half the slab. Above each profile 
panel we show the ratio of each model’s predictions to those of 
hml2_evoT_rt (i.e. dex of difference in the bottom panel). 

By comparing models grey_fixT_B and hml2JixT_B we 
can isolate the effect of using the grey approximation. The 
monochromatic model is characterized by a fully neutral hydro¬ 
gen core of approximately 1.5 kpc and sharp edges. Using the 
polychromatic spectrum smooths the hydrogen ionization pro¬ 
file and reduces the neutral hydrogen column density through 


7. Stromgren Sphere Solver 

rabacus has two sphere classes: SphereStromgren which 
models a point source at the center of a sphere and Sphere- 
Bgnd which models a sphere placed in a uniform and isotropic 
background. There is no closed form solution that describes the 
ionization profile in either case, however solutions are readily 
found using numerical integration. In this section and the next 
we present the solutions of rabacus and compare them to the 
results of the photo-ionization code cloudy version 13.03 last 
described in |Ferland et al.[ ( j2013| . 
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Figure 3: Ionization fraction and temperature profiles in five simple models of a Lyman limit system. Each panel shows the ionization fraction of one ion or the 
temperature. The five models considered are: 1) grey_fixT_B (thin solid) in which the grey approximation to the HM12 spectrum is used, the temperature is kept 
fixed, and recombinations are treated using case B rates, 2) hml2_fixT_B (thin dashed) in which the full HM12 spectrum is used, the temperature is kept fixed, 
and recombinations are treated using using case B rates, 3) hml2_evoT_B (thick dashed) in which the full HM12 spectrum is used, equilibrium temperatures are 
calculated, and recombinations are treated using case B rates, 4) hml2_evoT_A (dotted) as in hml2_evoT_B but with case A recombination rates (i.e. recombination 
radiation is ignored) and 5) hml2_evoT_rt (thick solid) in which the full HM12 spectrum is used, equilibrium temperatures are calculated, recombination rates are 
case A, and recombinations are treated by calculating the transfer of radiation through the system. Above each profile we show the ratio of each model’s prediction 
to those of hml2_evoT_rt (i.e. dex of difference in the profile panels). The scale in the top panels is a factor of ten in either direction. 



























We do not expect perfect agreement between the two codes 
for several reasons. First, cloudy makes use of a more detailed 
physical model including multi-level atoms and a large number 
of line transfer processes which are not included in rabacus. 
Second, the convergence criteria for cloudy allow for a devia¬ 
tion from unity in the sum of the ionization fractions of a given 
element (for example x m + v HII ) at the level of 10 -4 . Third, 
the atomic rates that determine the exact ionization and tem¬ 
perature structure in a parcel of gas are only known to a few 
significant digits at best and are not identical in the two codes. 
This is why we have tested rabacus in previous sections against 
closed form analytic solutions. It follows that differences be¬ 
tween cloudy and rabacus should be considered a measure of 
uncertainty as opposed to an error in either code, both of which 
have been verified against analytic solutions. 

The transfer of radiation from a central point source (ignor¬ 
ing recombination radiation for the moment) is a one-dimensional 
problem. This means the same techniques that we used in the 
case of plane parallel radiation incident on a slab geometry can 
be applied to the sphere as long as we include the inverse square 
geometric dilution of the radiation. In the following sections we 
verify the spherical solver in rabacus using a simple test prob¬ 
lem and then study the effects of temperature evolution, helium, 
and recombination radiation. 


7.7. Pure Hydrogen, Fixed Temperature 

We begin by examining Test 1 described in Iliev et al. ( [2006 ) 
which considers a pure hydrogen sphere of radius 6.6 kpc with 
uniform density, n u = 10 -3 cm -3 , and temperature, T = 10 4 K. 
The source has a photon luminosity of L n = 5.0 x 10 48 per sec¬ 
ond. The test specifies a monochromatic spectrum with all pho¬ 
tons having energy 13.6 eV. However, all spectra in cloudy have 
a finite width and so creating a monochromatic spectrum at the 
Lyman limit produces some non-hydrogen ionizing photons. In 
order to make a clean comparison to rabacus we instead use a 
monochromatic spectrum at 1.1 Ry or 15 eV. Recombination 
photons are treated by using the case B recombination rate (i.e. 
the on-the-spot approximation). In Figure [4] we present ioniza¬ 
tion profiles from both codes. The relative difference between 
the two solutions is typically less than one percent but can be a 
few percent near the ionization front. 


7.2. Temperature Evolution 

The second test from Iliev et al] ( [2006] ) allows for temper¬ 
ature evolution and includes a black body spectrum with effec¬ 
tive temperature T e ff = 10 5 K. We only consider photons with 
energy between 1 and 10 Rydbergs and normalize the spectrum 
such that the source emits L n = 5 x 10 48 photons per second. 
The temperature profile of the sphere is determined by balanc¬ 
ing photo-heating and radiative cooling processes. We present 
the solutions from rabacus and cloudy in Fig. [5] The relative 
difference in ionization fractions and temperature between the 
two codes for this test is on the order of a few percent except in 
the innermost quarter kpc. These differences coincide with the 
radius at which x m approaches the magnitude of the deviation 
from unity of the ionization fractions, Su = \l-x } 
cloudy model. 


x HII |, in the 




<1 






Figure 4: Ionization fractions for Test 1 of Iliev et al. ( 2006) (with 15 eV pho¬ 
tons). In the lower panel we show ionization fractions from rabacus (thick 
lines) and cloudy (thin lines). The deviation of the sum of the ionization frac¬ 
tions from unity (dn = 11 - x HI - x HII |) in the cloudy model is shown as a thin 
grey line. In the middle panel we show the absolute difference in ionization 
fractions cloudy - rabacus. In the top panel we show the relative difference in 
ionization fractions ( cloudy - rabacus ) / rabacus. In the bottom panel, the 
rabacus solution is shown in black and the cloudy solution is shown in green. 
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7.3. Inclusion of Helium 





Figure 5: Ionization and temperature profiles for Test 2 of Ilievet al~] (2006) . In 
the upper three panels we show ionization fractions from rabacus (thick black 
lines) and cloudy (thin green lines). The deviation of the sum of the ionization 
fractions from unity (<5 h = |1 - x m - x HII |) in the cloudy model is shown as a 
thin grey line. The top two panels show the absolute difference in ionization 
fractions, cloudy - rabacus, and the relative difference in ionization fractions, 
( cloudy - rabacus ) / rabacus. The bottom three panels show the same thing 
for temperature. The difference between the two codes in the temperature pro¬ 
file panel is only significant for the innermost quarter kpc. These differences 
coincide with the radius at which x HI approaches 6u in the cloudy model. 


Friedrich et al. ( 2012j ) extend the spherically symmetric tests 
described in Iliev et al.[ |2006| ) to include a uniform helium den¬ 
sity of n He = 8.7 x 10 -5 cm -3 and increase the radius of the 
sphere to 20 kpc. To verify the treatment of helium in rabacus 
we use a 10 5 K blackbody spectrum, case A recombination rates 
(i.e. ignore recombination radiation), and allow temperatures to 
vary as in the previous test. We consider photons with energy 
between 1 and 10 Rydbergs and normalize the spectrum such 
that the source emits L n = 5 x 10 48 photons per second. Except 
for the addition of temperature evolution, this is the same setup 
as Test 1 A in [Friedrich et al.| ( |2012| ). We use both rabacus and 
cloudy to generate solutions and present the ionization and tem¬ 
perature profiles in Fig. [6] At small radii we see the same effect 
as in the previous test due to the neutral fractions approaching 
S n = 11 - x HI - v HII | and d He = |l-x HeI -A HeII -A HeIII |. At larger radii, 
the relative differences in the hydrogen and helium neutral frac¬ 
tions are always less than 5%. In the case of the ionized species, 
the relative differences between ionization fractions in rabacus 
and cloudy is less than 10% for most radii. However, for ion¬ 
ized hydrogen and doubly ionized helium, the differences can 
be greater. It is hard to determine the exact cause of these dif¬ 
ferences, but the most likely explanations are differences in the 
atomic rates used in the two codes and line processes in cloudy 
that are not included in rabacus. The absolute difference be¬ 
tween ionization fractions in the two codes is always less than 
0 . 02 . 


7.4. Hydrogen Recombination Radiation 

In the previous Stromgren sphere examples, we have ei¬ 
ther ignored recombination radiation (i.e. used case A rates) 
or treated it using the on-the-spot approximation (i.e. used case 
B rates). A more accurate treatment would allow photons pro¬ 
duced via recombinations to travel some distance in the system 
before producing a photo-ionization. The outward only approx¬ 
imation is commonly used in analytic treatments of Stromgren 
spheres (e.g. [Ritzerveld 2005). In this approximation, all re¬ 
combination photons are assumed to be emitted in the outward 
radial direction. This is equivalent to a situation in which half 
the photons are emitted radially inward and half radially out¬ 
ward if the optical depth seen by recombination photons inte¬ 
rior to the emitting layer is zero. |Raicevic e t al. ( 2014) examine 
the effects of isotropically transporting recombination radiation 
in the Stromgren sphere problem using the same physical setup 
as in Test 1 of Iliev et al. ( |2006[ ). rabacus includes options to 
treat the isotropic emission of recombination radiation as well 
as options to use the outward only approximation. In this sec¬ 
tion, we present solutions produced by rabacus with isotropic 
transport of recombination radiation. 

When radiation is not constrained to move in the radial di¬ 
rection, one must account for the angular dependence of the op¬ 
tical depth when calculating photo-ionization and photo-heating 
rates. In spherical geometries this angular dependence is re¬ 
stricted to a dependence on the polar angle /a = cos 6. To calcu¬ 
late the transfer of recombination radiation, one must calculate 
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Figure 7: Illustration of ray tracing in spherical geometries. Because of sym¬ 
metry, all rays (dashed line) can be rotated such that they are perpendicular to 
the vertical axis. We indicate the origin of the rays with a dot in the center of 
the / = 1 layer and the radial direction with a grey solid line. In the calculation 
of path length ds through each layer we label the edges (from large to small 
radii) with the integer e , the layers (from large to small radii) with the integer / 
and the ray segments (from left to right) with the integer j. 



0 5 10 15 20 
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Figure 6: Ionization and temperature profiles for Test 1 A of |Friedrich et al.| 
|2012) with the addition of temperature evolution. In the upper three panels 
we show ionization fractions from rabacus (thick lines) and cloudy (thin lines). 


the optical depth through the sphere along directions for ev¬ 
ery layer. Each ray begins at the center of a layer and the di¬ 
rections, Hi are chosen using Gauss-Legendre quadrature. The 
number of directions can be specified by the user but defaults to 
Np = 32. We found that increasing or decreasing the number of 
angles by a factor of 2 had no effect on the tests presented here, 
however users should check convergence for each problem. 

With the proper rotation, all rays can be made perpendicular 
to the vertical axis. In Fig. [7] we show two rays that originate 
from the layer 1=1. The first (extending to the left) represents 
a direction with Hi = cos 6 > 0 while the second (extending 
to the right) represents a direction with Hi = cos 6 < 0. We 
label the edges of the sphere (from large to small radii) with 
the integer e , the layers of the sphere (from large to small radii) 
with the integer /, and the segments of the ray (from left to 
right) with the integer j. Given a layer of origin and a direction 
Hi, we first calculate the impact parameter b and the layer of 
closest approach 4 (e.g. 4 = 2 in Fig. [7]). In addition we 
label the radius of each edge with the function r e (i). Given this 


The deviation of the sum of the ionization fractions from unity (bn = |1 - x m - 
x mi \ and dn e = |1 - v HeI - x HeII - jt Hem |) in the cloudy model are shown as 
thin grey and green lines respectively. The top two panels show the absolute 
difference in ionization fractions, cloudy - rabacus, and the relative difference 
in ionization fractions, ( cloudy - rabacus ) / rabacus. The lines from the two 
codes are only distinguishable for the ionized species at large radii. The bottom 
three panels show the same thing for temperature. The difference between the 
two codes in the temperature profile panel is only visible at small radii. 


information, the path length ds(j) through segment j can be 
calculated as follows, 


1 

(N 

1 

+ 

yjfeU) ~ b2 

if j < 4, 

yjr 2 e(j) - b 2 


if j = h. 

Jr 2 (k + 1) - b 2 - 

yjr 2 (k)-b 2 

if j > 4. 
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Figure 8: Results for the test described in [Raicevic et al.|(2014) . In the upper 
panel we show the neutral hydrogen fraction resulting from using case A rates 
(ignoring recombination radiation), using case B rates (using the on-the-spot 
approximation) and from transporting recombination radiation isotropically. In 
the lower panel we show photo-ionization rates due to the central source (thick 
lines) and due to recombination radiation (thin lines) all with geometric dilution 
scaled out. In the middle panel we show the ratio of diffuse to central source 
rates. 


When modeling an isotropic background external to the sphere, 
one adds the optical depth from each ray segment j between 
the ray origin and the surface of the sphere. This is then used to 
attenuate the incoming radiation. When modeling recombina¬ 
tion radiation, emission from each layer is attenuated using the 
appropriate optical depth segments. 

In the upper panel of Fig. [8] we show hydrogen ionization 
profiles when using case A recombination rates (i.e. ignor¬ 
ing recombination radiation), when using case B recombination 
rates (i.e. the on-the-spot approximation), and when transport¬ 
ing recombination radiation isotropically. As expected, the case 
A result produces the smallest ionized region. The case B re¬ 
sult assumes that each recombination produces a local ioniza¬ 
tion and so produces a larger ionized region. When recombi¬ 
nation photons are emitted isotropically, they travel through the 
low optical depths at small radii without being absorbed and 
the solution tracks the case A model. However, all photons are 
eventually absorbed in the Stromgren sphere problem and so 
the final radius of the Stromgren sphere must be the same as 
in the case B approximation due to photon conservation. This 
result does not hold for slab geometries which are optically thin 
at their center (see §6.31 ). 

An instance of the SphereStromgren class has attributes 
that store photo-ionization rates due to both the central source 
and to diffuse recombination radiation. In the lower panels of 
Fig. [8] we show how these scale with radius. The diffuse field 
for the on-the-spot approximation is calculated using what we 
term the effective photo-ionization rate. In addition, to track¬ 
ing the source and diffuse photo-ionization rates, rabacus will 
define an effective photo-ionization rate, T eff , for each ionized 
species if recombination photons are not being transported. The 
effective rates are defined for each layer to be those that would 
result in the same ionization fractions if case A recombination 
rates were used. By definition the effective photo-ionization 
rate is equal to the actual photo-ionization rate if case A rates 
are used. However, if case B rates are used, subtracting the 
actual photo-ionization rate from the effective photo-ionization 
rate leaves the diffuse photo-ionization rate. The on-the-spot 
approximation produces a constant ratio of approximately 0.6 
between diffuse and central source photons. When ray tracing 
is used to follow the isotropic emission of recombination pho¬ 
tons, their contribution peaks near a radius of 4 kpc on the flat 
part of the ionization profile and goes to zero at the center of 
the sphere. These results are consistent with the findings of 
[Raicevic et al.[|2014l ). 


7.5. Helium Recombination Radiation 

Few papers have focused on the effects of diffuse radiation 


in the presence of hydrogen and helium (although see Can- 


|talupo and PorcianT||2011| ). In rabacus, we allow for the trans¬ 
port of radiation due to recombinations to the ground state for 
all ions. This allows for meaningful comparison to case A and 
case B approximations. In future versions of the code we will 
include ionizing photons produced by other processes such as 
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Baimer series recombinations, two photon processes, and sec¬ 
ondary ionizations and heating from fast (E > 100 eV) elec¬ 


trons ((Shull and van Steenberg| |1985[ (Furlanetto and Stoever 

20TO). 

To examine the effects of helium recombination photons we 
create a point source with a photon luminosity of L n = 10 48 
photons per second and a powerlaw spectrum with L v oc y~ a 
between 1 and 60 Rydbergs. Fig. [9] is analogous to Fig. [8] but 
for helium ionization profiles and photo-ionization rates. For 
isotropic transport of recombination photons we find that the 
diffuse photo-ionization rate for neutral helium is never greater 
than 0.2 of the direct source rate but peaks around 4 kpc for 
singly ionized helium as in the case of hydrogen. While the 
case B approximation is more accurate for the outer radii in 
our previous tests involving hydrogen, the case A approxima¬ 
tion is more accurate for helium. This result is due to the fact 
that helium recombination photons are less effective at ionizing 
helium than hydrogen recombination photons are at ionizing 
hydrogen. In fact, the opacity seen by helium recombination 
photons comes mostly from neutral hydrogen. In addition, the 
higher energy and lower photo-ionization cross-section of he¬ 
lium recombination photons allow some to escape the system 
all together (see jj | Appendix A| ). However, we show in Fig. 10 
that these results are sensitive to the steepness of the central 
source spectrum. 

To illustrate the previous point, we calculate the ionization 
structure in the sphere using five powerlaw spectra with spectral 
slopes a evenly spaced between 0 and 4 all normalized to emit 
L n = 10 48 photons per second. We present the results in Fig. 
10 In general, steeper spectral slopes produce a larger ratio 
of low-energy to high-energy photons leading to an increased 
singly-ionized helium fraction at r < 4 kpc. This in turn leads 
to a higher rate of recombination for singly ionized helium and 
the increased importance of recombination radiation relative to 
central source photons for ionizing neutral helium in the in¬ 
ner region. As the spectrum flattens, diffuse radiation becomes 
decreasingly important for both neutral and singly ionized he¬ 
lium, however the effect is stronger for neutral helium. Fur¬ 
thermore, a flatter spectrum contains more high-energy pho¬ 
tons than steeper spectra, which then doubly-ionize helium at 
a higher rate. However, the optical depth of its recombination 


radiation is small (see ^Appendix A| ) and most of the diffuse 
radiation escapes from the system. 


7.6. Isolating Physical Processes 

As a final examination of the capabilities of the Sphere- 
Stromgren class we solve five models similar to those described 
in the slab geometry section (see §6.3| ) to gain some insight into 
the effects of different physical processes. The sphere we con¬ 
sider has a radius of 20 kpc and hydrogen and helium number 
densities of n u = 10 -3 cm -3 and n Ue = 8.7 x 10 -5 cm -3 . We 
consider a 10 5 K blackbody source and its grey approximation 
(17.8 eV) both normalized to emit 5 x 10 48 photons per second. 
The polychromatic models are prefixed with tle5 as opposed to 
hml2 to indicate this change in spectrum. 

In Fig. [IT] we show ionization and temperature profiles for 
all five models. Above each profile panel we show the ra- 




Figure 9: Results for a version of the test described in |Raicevic et al.H2014) 
extended to include helium. In the upper panel we show the ionization fractions 
for all three helium ions resulting from using case A rates (ignoring recombi¬ 
nation radiation), using case B rates (using the on-the-spot approximation) and 
from transporting recombination radiation isotropically. In the lower panels we 
show the ratio of diffuse to central source photo-ionization rates. 
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Figure 11: Ionization fraction and temperature profiles in five simple Stromgren sphere models. Each panel shows the ionization fraction of one ion or the 
temperature. The five models considered are: 1) grey_fixT_caseB (thin solid) in which the grey approximation to the 10 5 K blackbody spectrum is used, the 
temperature is kept fixed, and recombination rates are case B, 2) tle5_fixT_caseB (thin dashed) in which the full 10 5 K blackbody spectrum is used, the temperature 
is kept fixed, and recombination rates are case B, 3) tle5_evoT_caseB (thick dashed) in which the full 10 5 K blackbody spectrum is used, equilibrium temperatures 
are calculated, and recombination rates are case B, 4) tle5_evoT_A (dotted) as in tle5_evoT_B but with case A recombination rates (i.e. recombination radiation is 
ignored) and 5) tle5_evoT_rt (thick solid) in which the full 10 5 K blackbody spectrum is used, equilibrium temperatures are calculated, recombination rates are case 
A, and recombination radiation is transferred through the system as opposed to being absorbed on the spot. Above each profile we show the ratio of each model’s 
prediction to those of tle5_evoT_rt (i.e. dex of difference in the profile panels). The scale in the top panels is a factor of three in either direction. 
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Figure 10: The effect of the powerlaw index, a, on the helium ionization struc¬ 
ture. In the upper three panels we show the ionization fractions for neutral, 
singly ionized, and doubly ionized helium. We show the results for five power- 
law spectra, L v oc y~ a , with a ranging between 4 and 0 (in black, red, orange, 
lime, and grey respectively). In each of these panels, the radii for which the 
ionization fraction is above 0.5 are indicated with a filled region. In the lower 
panels we show the ratio of diffuse to central source photo-ionization rates. As 
the spectrum becomes steeper, the contribution from recombination radiation 
becomes more important. 


tio of each model’s predictions to those of tle5_evoT_rt (i.e. 
dex of difference in the bottom panel). By comparing mod¬ 
els grey_fixT_B and tle5_fixT_B we can isolate the effect of 
using the grey approximation. The monochromatic model is 
characterized by a steep transition from ionized to neutral hy¬ 
drogen at approximately 7 kpc which can be seen clearly in the 
jCjjh panel. As in the case of the slab, using the polychromatic 
thermal spectrum smooths the hydrogen ionization profile. We 
note that the two models reach the same x HII value at large radii 
where collisional ionizations dominate over photo-ionizations 
for hydrogen. The harder photons that escape the ionized hy¬ 
drogen region have a larger cross-section for absorption due to 
helium than from hydrogen. 

By comparing models tle5JixT_B and tle5_evoT_B we 
can isolate the effect of self-consistently calculating the tem¬ 
perature structure. The equilibrium temperature is greater than 
10 4 K in the region in which hydrogen is ionized, but falls be¬ 
low that temperature at larger radii. For hydrogen, this produces 
lower neutral fractions at small radii and lower ionized fractions 
at large radii. This is because the ionization state of hydrogen 
at radii r > 7 kpc is determined by collisional ionizations. For 
helium, the higher temperatures at small radii also reduce x HeI , 
but helium ionizing photons still reach radii where the tempera¬ 
ture is cooler than 10 4 K and the reduced absorption of photons 
at small radii due to the higher temperatures leads to larger ion¬ 
ized fractions for both x HeII and v HeIII in model tle5_evoT_B. 

The larger case A recombination rate in model tle5_evoT_A 
produces more neutral species for both hydrogen and helium. In 
model tle5_evoT_rt we calculate the transfer of recombination 
photons isotropically emitted from each layer. As in the case of 
the slab, the rt and case A solutions agree in the optically thin 
center of the sphere. Interestingly, the case B approximation 
is more accurate for hydrogen while the case A approximation 
is more accurate for helium. Overall, the variations we have 
examined lead to differences of approximately a factor of 2 for 
the hydrogen ionization profiles and a factor of 50% for the 
helium ionization profiles. 


8. Background Sphere Solver 


In this section we describe the SphereBgnd class which 
models a spherically symmetric gas distribution placed in an 
isotropic background radiation field. Even in the simplified 
case of monochromatic radiation incident on a uniform den¬ 
sity and temperature sphere, there are no analytic solutions for 
the ionization profile. In addition, the code we have used for 
comparison thus far, cloudy, does not support this geometry. 

This geometry has been studied in the context of mini-halo 
evaporation (e.g. | Shapiro et al.| |2004| a nd in the context of 
hydrogen absorption systems (e.g. Zheng and Miralda-Escude! 


[20021 ) however a direct comparison with either work would not 
be fruitful. In the first case, ionization profiles are not pre¬ 
sented. In the second it would be difficult to know if differences 
arise from the set of included physical processes, the method of 
solution, or the implementation in code. Our goal for this final 
section is not a validation of the SphereBgnd class but rather a 
description of results in the fight of the other validated pieces of 
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the rabacus. A complete validation of this geometry must await 
a dedicated comparison project. 

In what follows we will describe a simple model for a gaseous 
halo immersed in the cosmic UV background radiation. Dark 
matter halos have been shown to have universal density profiles 


(Navarro et al. 1997) characterized by a scale radius R s . To 
model the gaseous component of the halo we create an NFW 
profile with a core at r - 37? s /4 as in [Mailer and Bullock] ( |2004| ) 
and [Barnes and Haehnelt ( 2014 ). Specifically, the virial radius 
R y is related to the virial mass M v as, 


R y = 46.1 kpc 


Ay£l m h 2 

24.4 


-1/3 


1 + Z 
~ 33 ~ 


-l 


My 


10 n M o 


1/3 


(41) 


where A v is the mean overdensity inside the virial radius of 
the halo. The ratio of the virial radius and the scale radius is 
known as the concentration parameter c = Ry/R s . For our sim¬ 
ple model we take the average relationship between concentra¬ 
tion and virial mass described in [Macci6 et al.| ( |2007| , 
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with cq = 3.5. The gaseous hydrogen density profile is, 


n„(r) 


_ R&o _ 

[ r + (3/4)/? s ](r + R s ) 2 


(43) 


where po is a constant set such that the integrated hydrogen gas 
mass is equal to M H = fuMy and /h = (1 - T p )Qb/^m- The 
profile is then fully specified by choosing a virial mass, a mean 
overdensity inside the viral radius, and a redshift which we take 
to be My = 10 10 M o , A v = 200, and z = 0 for our example. This 
results in a concentration parameter of 18. 

The ray tracing tools necessary to calculate the transfer of 
the isotropic background through the spherical geometry are 
the same we described in the context of isotropic recombina¬ 
tion radiation. In other words, any ray we need to trace for 
treating an isotropic background can also, through the appro¬ 
priate rotation, be made to be perpendicular to the vertical axis. 
In Fig. [12] we show ionization and temperature profiles using 
the same five models as in Fig. [3j grey_fixT_B, hml2JixT_B, 
hml2_evoT_B, hml2_evoT_A, and hml2_evoT_rt. In this case, 
the steep density profile serves to lessen the differences between 
the models relative to the constant density slab case. This is 
partly due to the fact that the chosen halo mass produces a neu¬ 
tral core for all models. Models are most sensitive to changes 
when the core is near the optically thin to optically thick tran¬ 
sition. Examining in detail the contribution to absorber abun¬ 
dances from different halo masses is beyond the scope of this 
paper, but one of the primary intended applications of rabacus. 


9. Conclusion 

We have described rabacus, a Python package for calculat¬ 
ing the transfer of hydrogen ionizing radiation in simplified ge¬ 
ometries relevant to astronomy and cosmology. We presented 
example solutions for: 1) a semi-infinite slab gas distribution 


with isotropic radiation incident from both sides, 2) a spheri¬ 
cally symmetric gas distribution with a point source at the cen¬ 
ter, and 3) a spherically symmetric gas distribution in a uni¬ 
form background. We have also demonstrated the effects of 
temperature evolution, the inclusion of helium, and the treat¬ 
ment of recombination radiation on the ionization and tempera¬ 
ture profiles in these geometries. For slabs exposed to isotropic 
backgrounds at Lyman limit system densities, these processes 
can change hydrogen and helium column densities by an order 
of magnitude. In the case of spherical geometries the magni¬ 
tude of these effects is not as large but can still change ion¬ 
ization fractions by factors of 2. The flexibility and speed of 
rabacus lend itself to rapid prototyping of ideas and allows 
it to serve as a check for more complicated radiative transfer 
schemes. It is made publicly available and is written in an in¬ 
terpreted language, Python, with a large user base. Future im¬ 
provements will include the addition of metal line cooling and 
non-equilibrium solvers. 
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Appendix A. Recombination Radiation 

In this section we examine in more detail the Stromgren 
sphere test described in §7.5| Our goal is to determine the total 
opacity encountered by recombination photons and the relative 
contribution to the total from different ions. We will make use 
of the quantity r r which is simply the optical depth through 
each thin radial shell of thickness dr, 


T n.J dr = n m a ul (E m ) 

T\ A Jdr = n m cr m (E m ) + n HeI o- HeI (£ HeI ) 

t 54 J dr = "ffiO’ra^Hen) + + ^Hell ^Hell (£„en) (A.1) 


In Fig. A.13| we show the radial dependence of several quan¬ 
tities. The top panel shows the emission coefficients as calcu¬ 
lated in Eqs. [37]and[38] The majority of recombination photons 
emitted at any radii come from hydrogen. At inner radii, re¬ 
combinations of doubly ionized helium produce the majority of 
helium recombination photons while at radii greater than 5 kpc 
both types of helium recombination produce a similar number 
of photons. 

In the middle three panels we show how the total opacity 
is partitioned amongst the absorbing species. Photons at 13.6 
eV are only absorbed by hydrogen and so the optical depth pro¬ 
file follows that of the n m profile (see bottom panel in the same 
figure). Photons at 24.6 eV are absorbed mostly by neutral hy¬ 
drogen. The cross-section for absorption by helium is six times 
larger than that of hydrogen at this energy, but the greater den¬ 
sity of neutral hydrogen at all radii more than compensates for 
this. Photons at 54.4 eV can be absorbed by all three absorb¬ 
ing species. The cross-sections at this energy for both helium 


16 


























Figure 12: Ionization fraction and temperature profiles in five simple models of a gaseous halo in the UV background. Each panel shows the ionization fraction 
of one ion or the temperature. The five models considered are: 1) grey_fixT_caseB (thin solid) in which the grey approximation to the HM12 spectrum is used, the 
temperature is kept fixed, and recombination rates are case B, 2) hml2_fixT_caseB (thin dashed) in which the full HM12 spectrum is used, the temperature is kept 
fixed, and recombination rates are case B, 3) hml2_evoT_caseB (thick dashed) in which the full HM12 spectrum is used, equilibrium temperatures are calculated, 
and recombination rates are case B, 4) hml2_evoT_A (dotted) as in hml2_evoT_B but with case A recombination rates (i.e. recombination radiation is ignored), and 
5) hml2_evoT_rt (thick solid) in which the full HM12 spectrum is used, equilibrium temperatures are calculated, recombination rates are case A, and recombination 
radiation is transferred through the system as opposed to being absorbed on the spot. Above each profile we show the ratio of each model’s prediction to those of 
hml2_evoT_rt (i.e. dex of difference in the profile panels). The scale in the top panels is a factor of three in either direction. 
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species differ by less than ten percent, but singly ionized helium 
has a cross-section that is 13 times larger than that of hydrogen. 
At small radii, singly ionized helium absorbs the majority of 
photons. At larger radii, neutral hydrogen and helium absorb 
a similar amount of photons but the absorption by hydrogen 
is always equal to or greater than that by neutral helium. In 
Fig. |A.14| we examine changes to the hydrogen ionization pro¬ 
file caused by the inclusion of helium in §7.5| There are two 
competing effects to consider. Helium will absorb some of the 
source photons that would have been absorbed by hydrogen, 
however some fraction of the helium recombination photons 
will be absorbed by hydrogen. In the bottom two panels of Fig. 
|A.14| we show that the first is larger and that the neutral hy¬ 
drogen fraction increases by as much as twenty percent when 
helium is included. Part of the reason that the helium photons 
do not make a larger contribution is because they escape the 
system. This can be seen in the top panel in which we show 



(A.2) 



The optical depth to the surface of the sphere for photons at 
54.4 eV is only larger than unity at radii smaller than 6 kpc. 
While recombination photons at 24.6 eV do encounter signif¬ 
icant optical depth, the abundance of these photons is always 
lower than those from hydrogen recombinations (see top panel 
of Fig. [03}. 
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the Stromgren sphere test in ClU In the top panel we show r >r (see Eq. |A,2) 
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the bottom two panels we show the neutral hydrogen ionization fraction profile 
with and without helium. The inclusion of helium enhances x HI by as much as 
25%. 
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